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COMPARISON OF FINITE-DIFFERENCE SCHEMES 
FOR ANALYSIS OF SHELLS OF REVOLUTION 


By Ahmed K. Noor* and Wendell B. Stephens 
Langley Research Center 

SUMMARY 

Several finite- difference schemes are applied to the stress and free-vibration anal- 
ysis of homogeneous isotropic and layered orthotropic shells of revolution. The study is 
based on a form of the Sanders- Budiansky first-approximation linear shell theory modi- 
fied such that the effects of shear deformation and rotary inertia are included. A Fourier 
approach is used in which all the shell stress resultants and displacements are expanded 
in a Fourier series in the circumferential direction, and the governing equations reduce 
to ordinary differential equations in the meridional direction. 

While primary attention is given to finite-difference schemes used in conjunction 
with first-order differential- equation formulation, comparison is made with finite- 
difference schemes used with other formulations. These finite-difference discretization 
models are compared with respect to simplicity of application, convergence character- 
istics, and computational efficiency. Numerical studies are presented for the effects of 
variations in shell geometry and lamination parameters on the accuracy and convergence 
of the solutions obtained by the different finite- difference schemes. 

On the basis of the present study it is shown that the mixed finite -difference scheme 
based on the first-order differential- equation formulation and two interlacing grids for the 
different fundamental unknowns combines a number of advantages over other finite- 
difference schemes previously reported in the literature. 


INTRODUCTION 

In the past decade a substantial capability has been developed for the numerical anal- 
ysis of rotationally symmetric shells. The most widely used numerical analysis proce- 
dures for this class of shells are those based on numerical integration, finite elements, 
and finite differences. A number of publications exist which review and assess the rela- 
tive merits of these three numerical techniques (see, for example, refs. 1, 2, and 3). 
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The large majority of existing finite-difference programs for rotationally sym- 
metric shells are based on either the displacement formulation (with the fundamental 
unknowns being displacement parameters) (ref. 4), or the Budiansky-Radkowski second- 
order differential-equation formulation (refs. 1, 5, and 6). While the advantages of using 
a first-order differential-equation formulation have been widely recognized in numerical 
integration techniques and in matrix progression method (refs. 7 to 11), only limited use 
of this formulation has been made in the case of finite differences (refs. 12 and 13). 

The objective of this paper is to assess the relative merits of several finite- 
difference schemes used for the linear elastic stress and free-vibration analysis of 
rotationally symmetric shells. Primary attention is given to mixed finite- difference 
schemes used in conjunction with the first-order ordinary-differential-equation formula- 
tion. However, some consideration is given to finite-difference schemes used with the 
displacement formulation and with the Budiansky-Radkowski second-order differential 
equations. 

The term "mixed" refers to the fact that both stress resultants and displacements 
are chosen as primary variables. The analytical formulation is based on the Sanders- 
Budiansky linear shell theory. The first-order ordinary-differential-equation formula- 
tion of that theory is modified such that the effects of shear deformation and rotary inertia 
are included. 
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SYMBOLS 

dimensionless elastic compliance coefficients of the shell 

dimensionless elastic coefficients of the shell 
dimensionless elastic stiffnesses of the shell 
reference extensional rigidity of the shell 
elastic stiffness of the kth layer of the shell 
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elastic modulus for isotropic materials 
error index (see eq. (19)) 

elastic moduli in the direction of the fibers and normal to it, 
respectively 

shear moduli in the plane of the fibers and normal to it, respectively 
local thickness of the shell (see fig. 1) 

distance from the reference (middle) surface to the top and bottom 
surfaces of the kth layer, respectively 


L 

l 


m 


NL 

n 

Po 
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reference length of the shell 

finite- difference interval in the meridional direction 
moment stress resultants (see fig. 2) 

Fourier harmonic in the circumferential direction 
direct stress resultants (see fig. 2) 
modified (boundary) stress resultant 
symmetric stress resultants 

number of layers of laminated shells 

number of finite -difference intervals in the meridional direction 
vectors of external forces 
amplitude of load harmonic 

intensity of external loading in the coordinate directions (see fig. 2) 
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Q^,Qo 

R 

RpRg 


transverse shear stress resultants (see fig. 2) 
radius of toroidal cross section 

principal radii of curvature in the meridional and circumferential 
directions, respectively 

normal distance from the shell axis to the reference surface 
(see fig. 1) 


time 


x l’ x 2’ x 3 

k,K 1 


displacement components of the reference surface in the coordinate 
directions (see fig. 2) 

curvilinear coordinate system (see fig. 1) 

dimensionless frequencies 


•Poisson's ratio for isotropic materials 
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Poisson's ratio measuring strain in T-direction due to uniaxial normal 
stress in the L-direction 

dimensionless meridional and circumferential coordinates, 
respectively (see fig. 1) 

density of the kth layer of laminated shells 
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nondimensional geometric parameters 
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rotation components (see fig. 2) 


<Pg'<Pe 

{^m} anc * {^ m } vectors of fundamental unknowns defined by equations (14) and (15) 

a circular frequency of vibration of the shell 

LI row matrix 

{ } column matrix 

[ ] rectangular or square matrix 


Subscripts: 

i generic finite-difference station 

L,T denote the direction of fibers and the transverse direction 

m denotes the mth Fourier harmonic 

Finite-difference models: 

DW displacement formulation, whole- station scheme 

SW second-order formulation, whole- station scheme 

MW mixed formulation, whole- station scheme 

MHS 1, MHS 2 mixed formulation, half -station scheme, single grid 

MHI mixed formulation, half- station scheme, interlacing grid 

MATHEMATICAL FORMULATION 
Shell Geometry 

Figure 1 shows the geometric characteristics of a rotationally symmetric shell as 
follows: h is the local thickness of the shell; r is the normal distance from the shell 
axis to the reference surface; Rj and Rg are the principal radii of curvature in the 
meridional and circumferential directions, respectively. The expressions for Rj^ and 
Rg are given in reference 5. 
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To cast the problem in nondimensional form, the two dimensionless coordinates 
i and 6 are used, where 


£ 



(1) 


9 



( 2 ) 


in which x^ and Xg are the meridional and circumferential coordinates, respectively, 
and L is a reference length. 

Also, as in reference 5, the following dimensionless geometric parameters are 
introduced: 


P 



(3) 





(4) 


and 


L 



P 


(5) 


( 6 ) 


Reduction to a One-Dimensional Problem 

The analytical formulation is based on a form of the Sanders- Budiansky linear 
theory modified such that the effects of shear deformation and rotary inertia are included. 

If the shell stress resultants and displacements as well as the external loads are 
expanded in a Fourier series in the circumferential direction, then 

00 

[_ U # W N f M € «f] = E Km W m *i, rn N «,m M f,m %,m] cos < 7 > 

m=0 
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L U ^ % % Q*] = 2 L\m N ^,m M ^,m Q ,,m] sin mr (8) 

m=0 

00 

[_ p £ p r] = H h,m p r,m] cos m 0 
m=0 L 

and 

00 

p e = H p 8,m sin m s (10 > 

m=0 

where 



in which and are the modified (symmetric) stress resultants defined in 

reference 5. 

For free-vibration analysis, the right-hand sides of equations (7) and (8) must be 
multiplied by ^ where cc is the circular frequency of vibration of the shell. 


Governing Differential Equations 

The stress and free-vibration problems for shells of revolution can be formulated 
in a number of different ways. For example, such problems can be formulated in terms 
of (1) the three displacement components of the reference surface, (2) coupled second- 
order ordinary differential equations in the generalized displacements and meridional 
moment stress resultant, or (3) first-order ordinary differential equations in the shell 
stress resultants and displacements. In the absence of shear deformation, the governing 
equations of the displacement formulation for cylindrical shells are given in reference 14, 
and the governing equations of the second-order formulation for general shells of revolu- 
tion are presented in references 1 and 5. These equations are not reproduced here. On 
the other hand, the first-order equation formulation which has been extensively studied 
herein is outlined below. Discussion is focused on the choice of the fundamental unknowns 
and the form of governing differential equations in such a way as to enhance the compu- 
tational efficiency of the finite-difference discretization. 

The fundamental unknowns of the first-order differential-equation formulation used 
herein are chosen to be the shell quantities that appear in the statement of the boundary 
(or interface) conditions along a parallel circle. These include: the five generalized 
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displacements U^, U 0 , W, <pg, and <p Q -, and the five stress resultants N^, N^, 

Q^, ML, and (see fig. 2). 

By properly ordering the fundamental unknowns, the governing differential equa- 
tions for the mth Fourier harmonic can be arranged to yield a symmetric coefficient 
matrix of the fundamental unknowns as follows: 
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The term C Q is a reference extensional rigidity of the shell. The shell stiffnesses and 
compliances are given in appendix A. The [s a ], [s b ],and [s] are 5 x 5 matrices, the 
first two matrices are symmetric ; superscript T denotes transposition; the matrix [Ijj 
is a diagonal matrix; [m a ] and [m^] are symmetric mass matrices; and {P^ } and 

{Pg m j are vectors of external forces. The formulas for the coefficients of the afore- 
mentioned matrices are given in appendix B. For stress -analysis problems co = 0 and 
for free -vibration problems {P^ = {Pg = 0. 

The choice of the fundamental unknowns shown in equations (14) and (15) is similar 
to that suggested by Goldberg (ref. 7); however, the present formulation utilizes the geo- 
metric parameter p (^) in defining the unknowns for stress resultants. This particular 
choice allows casting the governing equations in a symmetric form and minimizes the 
number of nonzero terms in the matrix [s]. 
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FINITE -DIFFERENCE DISCRETIZATION 


The finite-difference models used with the displacement and second-order 
differential-equation formulation were based on the whole-station approximation. These 
two models will be referred to subsequently as DW (displacement, whole -station) and 
SW (second-order, whole -station) schemes. The latter model, SW, was first suggested 
by Budiansky and Radkowski (ref. 5) and developed into an operational program "SALORS" 
in reference 15. Since the DW and SW schemes have been discussed in references 14 and 
15, the discussion is confined here to the discretization of the first-order differential 
equations. To this end, it is convenient to express equations (12) as follows: 

M { H m} + [#m} + Dllbm} = -“V] { H m} ' { P l,m} (16a) 

and 

[ sb ]{*m} + [ S ] T { H m}-[ I l]{0 H m} = -‘' 2 [ mb K*m}-{ P 2,m} < 16b > 

Four finite -difference models have been used for the discretization of equations (16). 
These models differ by the difference -quotient expressions used for approximating the 
first derivatives and by the location of the points at which the difference equations are 
applied. The four models can be identified as mixed whole station (MW), mixed half- 
station, single grid (MHS 1 and MHS 2), and mixed half- station, interlacing grid scheme 
(MHI). The two schemes MHS 1 and MHS 2 differ one from the other in the difference- 
quotient expressions used for approximating the first derivatives. For shells with con- 
stant elastic and geometric characteristics, the discrete models obtained by these two 
schemes are identical. 

The characteristics of the four finite-difference models MW, MHS 1, MHS 2, and 
MHI along with those of the two models DW and SW are summarized in table 1. While 
three of the four schemes used with the first-order equation formulation, schemes MW, 
MHS 1, and MHS 2, have been reported previously in the literature, the fourth scheme 
(MHI) which was developed by the authors represents the prime thrust of the present 
work and is described in detail below. It should be noted, however, that other inter- 
lacing grids have been suggested and used previously by others (e.g., refs. 4 and 16). 

The basic idea of the mixed half-station interlacing-grid scheme (MHI) used in the 
present study is to define the first derivatives of each of the fundamental unknowns at 
points lying midway between the points of definition of the same unknowns. This can be 
accomplished by using two sets of interlacing grids for the two groups of fundamental 
unknowns |H m j and The quantities jV m | and |dH m j in equations (16) are 
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defined at points lying midway between the points of definition of {H m j and {d^ m }. 

For convenience both {H m | and {0 m } are defined at the boundaries and interfaces. 

The finite -difference equations which simulate the governing equations (eqs. (16)) 
are obtained by replacing the first derivatives in these equations by an appropriate 
difference-quotient expression, depending on the order of approximation desired (see 
ref. 17). 

Also, due to the coupling between and in equations (16a) and be- 

tween {H m } and {dH m j. in equations (16b) (i.e., due to the presence of the matrix 
[S] ), the values of {'Z' m } in the first equation and {H m j in the second equation have 
to be obtained through interpolation between their values at their respective control 
points. As will be shown subsequently, the accuracy of solutions is less sensitive to this 
interpolation than to the averaging of the difference-quotient expressions, which is used 
in conventional schemes (with all the fundamental unknowns and their derivatives defined 
at the same set of points). 

In order to maintain the number of unknowns equal to the number of equations, only 
the five equations (16a) are applied at the boundary, the other five equations are replaced 
by the five boundary conditions. The accuracy of the modified scheme presented herein, 
however, was found to be insensitive to the particular choice of the five equations to be 
applied at the boundary. A summary of the equations applied at the different nodal points 
is given in figure 3. Central differences are used at interior points and backward (or for- 
ward) unsymmetric differences are used at points lying on or near the boundaries (or 
interfaces). The number of such specialized unsymmetric difference formulas is de- 
pendent on the order of the interior discretization error used in the finite-difference 
scheme (see ref. 17). 

As an illustration, the difference equations are shown for the MHI scheme. These 

O 

difference equations of order 0{t ), where l is the finite -difference interval, at generic 
interior points i and i+1, are given by 


Kb * [t m 4 *4 pj * ; MilMw 


2 

2 


Mi { H m}i - { P l,ml i (17a) 


t sb ] i+ i {v)i + l + It Pll + \ { H m}i + [t + 1 MlJ { H m}i + 2 

= - 2 M i+ l{Mi + l - { p 2 ,m} i+1 < 17b > 


10 


The resulting finite-difference field equations can be represented in the following 
compact form: 

[ K m]{ Z m} = -" 2 [“H Z m}-{Pm} < l8 > 

where [K m ] and [rft] contain the generalized stiffness and mass distribution, |z m | 
is the vector of unknowns composed of the subvectors { H m }i and at the 

various finite -difference stations and |P m ] is the vector of external forces and thermal 
effects. Equations (18) are banded and their bandwidth depends on the orders of the in- 
terior and boundary approximations. The use of interlacing grids results in reducing 
both the interior discretization error and the bandwidth of the resulting finite-difference 
equations and, therefore, the computational efficiency of the scheme is improved. That 
this is indeed so can be seen from figure 4, where a schematic representation is given of 
the finite-difference equations corresponding to the four schemes MW, MHS1, MHS2, and 
MHI. 

The four finite -difference schemes MW, MHS 1, MHS 2, and MHI have a number of 
advantages also shared with other numerical approximation techniques based on the first- 
order equation formulation. These include the simplicity of the form of the governing 
differential equations, the absence of the derivatives of the elastic characteristics of the 
shell in these equations, and the simplicity of numerical discretization and of handling the 
boundary and interface conditions. As a result of the cited advantages, the effort required 
in the computer implementation of these difference schemes (coding, debugging, and veri- 
fication) has been found to be significantly reduced. 

EIGENVALUE EXTRACTION TECHNIQUE 

For free -vibration problems, a variant of the inverse-power method with shifts 
similar to that presented in reference 1 is used for the determination of the natural fre- 
quencies and mode shapes. In the present study, advantage was taken of the banded form 
of the matrix [kJ and a direct, Gaussian elimination procedure was used for each 
iteration to evaluate the new trial vector |z m j. 

NUMERICAL STUDIES 

In an attempt to assess the relative merits of the different finite -difference schemes, 
a large number of stress and free-vibration problems of homogeneous isotropic and layered 
orthotropic shells of revolution were solved using the aforementioned six finite-difference 
schemes. 
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The various finite-difference models are compared with respect to simplicity of 
application, convergence characteristics, and computational efficiency. The results of 
four typical problem sets representing different shell geometries and wall construction 
are discussed with respect to convergence characteristics herein. These problems are: 

a. Free vibrations of isotropic cylindrical shell 

b. Asymmetric stresses in isotropic toroidal shell subjected to a uniform normal 

pressure 

c. Asymmetric free vibrations of layered orthotropic toroidal shell 

d. Free vibrations of layered orthotropic spherical shells 

The four difference models based on the first-order equation formulation (MW, 

MHS 1, MHS 2, and MHI schemes) were applied to all four problems, and the second-order 
equation model (SW scheme) was applied to the first two problems. Comparison was made 
with the results obtained using the displacement model (DW scheme) for the first problem 
(ref. 14). 

In the first two problems in order to provide a meaningful comparison with finite- 
difference schemes SW and DW, the tracing constants k r and k a (where a = 1,2) in 
the appropriate equations in appendix B, were set equal to 0 and 10^, respectively. This 
resulted in dropping the rotary inertia terms and suppressing the shear deformation in 
these problems. 


Free Vibrations of Cylindrical Shells 

The first problem considered was that of asymmetric free vibrations of isotropic 
cylindrical shells with clamped edges. The characteristics of the shells considered are 
shown in figure 5. This problem was taken from reference 14 where it was concluded 
that the accuracy of the finite-difference method, based on a differential-equation formu- 
lation in terms of the three midsurface displacements, deteriorates rapidly as the length 
ratio L/R increases or the thickness ratio h/R decreases. Therefore, in order to 
assess the accuracy and rate of convergence for the six finite-difference schemes, the 
long thin shell shown in figure 5 was selected. For this problem it was found that both 
the shear deformation and rotary inertia were negligible and, therefore, the use of 
schemes DW and SW is justified. 

The nondimensional minimum frequencies obtained by the different finite-difference 
schemes are summarized in table 2 and the corresponding mode shapes and modal stresses 
obtained by schemes MHS 1, MHS 2, and MHI with 20 finite-difference intervals are shown 
in figure 6. 
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As is seen from this figure the three schemes MHS 1, MHS 2, and MHI accurately 
predict the mode shape and the modal stress-resultant N^, but for the modal moment 
M^, scheme MHI misses the spike at the boundaries and schemes MHS 1 and MHS 2 give 
an oscillatory type of variation. 

The results presented in table 2 show that the frequencies obtained by the mixed 
half-station interlacing grid finite -difference scheme (MHI) are more accurate than the 
corresponding frequencies obtained by other schemes (for the same number of finite- 
difference intervals). The error in the frequency obtained by using MHI scheme with 
10 finite-difference intervals is only 1.5 percent. On the other hand, the frequencies 
obtained by DW and SW schemes (based on the displacement and the four second-order 
equations) are considerably less accurate than those obtained by all other schemes based 
on the use of first-order differential-equation formulation. Also, the frequencies ob- 
tained by DW and SW schemes are identical to four significant digits. This is because 
the discretized equations in the SW scheme for cylindrical shells can be reduced to sets 
of equations identical to those obtained by direct discretization of the three differential 
equations of the displacement formulation, DW scheme (see ref. 17). Similar reasoning 
can be used to explain the fact that the frequencies obtained by MHS 1 and MHS 2 schemes 
are identical and those obtained by the MHI scheme are identical to the ones obtained by 
the MW scheme using twice as many finite -difference intervals (compare the correspond- 
ing entries in table 2). 


Stress Analysis of a Closed Toroidal Shell 

As a second example consider the asymmetric stress analysis of an isotropic closed 
toroidal shell subjected to a normal pressure which is uniform in the meridional direction. 
The characteristics of the shell and loading are shown in figure 7. 

Solutions have been obtained using the finite-difference schemes SW, MW, MHS 1, 
MHS 2, and MHI for the asymmetric stress analysis for an internal pressure loading which 
is uniform in the meridional direction. Only the second Fourier harmonic in the circum- 
ferential direction is considered (m * 2). For this shell and loading, the shear deformation 
was found to be negligible and therefore the use of scheme SW is justified. Due to sym- 
metry of the shell and loading, only half the meridian was considered and the symmetric 
boundary conditions at the ends are = 0. 

As a quantitative measure of the relative accuracy of the stress resultants and dis- 
placements obtained by the different finite -difference schemes, the following error index 
Ef (a function of f) has been introduced: 
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( 19 ) 



where 

f stands for any of the stress resultants or generalized displacements 

^ and f^ are the exact and approximate values, respectively, of the function at 
the ith finite-difference station. The exact value is taken to be the converged 
solution, and the approximate value refers to the value obtained by each of the 
finite-difference schemes. 

If I is the maximum absolute value of the exact (or converged) function in the 

l max l 

interval of interest (half the meridian of the shell in this case). 

n+1 is the total number of finite-difference stations used in the approximate 
solution. 

In order to simplify comparison with other schemes, the values of the unknowns 
{ i/j m } in scheme MHI were computed at the same nodal points as the unknowns [H m |* 
This was done after obtaining the modified solution. As an example, at the node point i 

{^m}i " ^ {{^mji-l + {^mji+l} 

where and are the values of jV m j at points i*l obtained from 

the modified solution. 

While the error index Ej (eq. (19)) is similar to the root-mean- square error, it 
has the added advantage that it gives less weight to the smaller values of the function, 
which are usually of less practical interest. The smaller the error index E^, the more 
accurate the approximate solution (obtained by the difference scheme) is. 

The values of the error index Ej for each of the stress resultants and generalized 
displacements obtained by the different finite -difference schemes using 10, 20, 30, and 
40 intervals in half the shell meridian are summarized in table 3. Also, figure 8 shows 
plots for the stress resultants and displacements obtained by schemes SW, MW, and MHI. 
Also shown are the essentially exact results (converged solution) for this problem. Solu- 
tions obtained by schemes MHS 1 and MHS 2 are not shown in order to avoid complicating 
the figure. To provide a reasonable comparison among the three finite -difference solu- 
tions shown in figure 8, the total number of algebraic equations for each method was made 
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comparable. This was achieved by using 10 finite-difference intervals in half the shell 
meridian in schemes MW and MHI and 20 intervals for scheme SW. Note that if the shear 
deformation is neglected, the governing differential equations used in schemes MW and 
MHI reduce to eight and the total number of algebraic equations for each method will 
almost be identical. Figure 8 and table 3 clearly show the high accuracy of the predic- 
tions of the proposed scheme (MHI) and that this accuracy occurs for both displacement 
and stress quantities in this case. 

Moreover, it was found that the accuracy of the solutions obtained by scheme MW 
was very sensitive to the particular choice of the five equations applied at the boundary. 
The best accuracy was obtained when the boundary conditions were applied in a manner 
similar to that used in finite-element method, wherein the 10 equations (eqs. (12)) are 
applied at the boundary, then the boundary conditions replace the equations with the 
specified boundary quantities (in this case displacement, rotation components, and stress 
resultants) appearing along the diagonal. The results shown in table 3 for scheme MW 
were obtained for the aforementioned choice. On the other hand, the accuracy of the 
solutions obtained by scheme MHI was found to be insensitive to the choice of the 5 equa- 
tions to be applied at the boundary. 


Asymmetric Free Vibrations of Laminated Orthotropic Toroidal Shell 

In order to study the effect of material orthotropy of the shell on the accuracy and 
rate of convergence of the different finite -difference schemes, the free -vibration problem 
of an eight-layered, graphite -epoxy composite shell was studied. The shell had the same 
dimensions and total thickness as that of the isotropic shell considered previously (see 
fig. 7). Both shear deformation and rotary inertia were accounted for and consequently 
scheme SW was not used. 

An indication of the accuracy and rate of convergence of the minimum frequency 
obtained by schemes MW, MHS 1, MHS 2, and MHI is shown in figure 9. The results ob- 
tained by schemes MHS 1 and MHS 2 were almost identical. As shown in figure 9 the 
frequencies obtained by all schemes converged to the same value, but those obtained by 
scheme MHI had a faster rate of convergence. 


Free Vibrations of Layered Composite Spherical Shell Segment 

As a final example consider the free vibrations of an eight layered, graphite -epoxy 
spherical-shell segment with clamped edges. The fibers of the different layers alternate 
between the circumferential and meridional directions, with the fibers of the top layer 
running in the meridional direction. Shells with different rise ratios, ranging from 
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shallow to deep shells, have been studied. The characteristics of these shells are shown 
in figure 10. For all these shells solutions were obtained using the first order difference 
models (schemes MW, MHS 1, MHS 2, and MHI). 

The solutions obtained using the different schemes are summarized in table 4. As 
can be seen from this table scheme MHI is superior to all other schemes. The accuracy 
of the solutions obtained by this scheme does not deteriorate as the shell becomes 
shallower. 


CONCLUDING REMARKS 

A comparison is made between a number of finite -difference schemes for analysis 
of shells of revolution. Primary attention is given to finite-difference schemes based on 
the use of first-order ordinary differential -equation formulation; however, some consid- 
eration is given to both the displacement and the second-order differential -equation models. 
The various finite-difference discretization models are compared with respect to simplic- 
ity of application, convergence characteristics, and computational efficiency. 

On the basis of the present study it is shown that the mixed finite -difference schemes 
have a number of major advantages in common with other numerical approximation tech- 
niques based on the first-order differential-equation formulation. These include the sim- 
plicity of the form of the governing differential equations, the absence of the derivatives 
of the elastic characteristics of the shell in these equations, and the simplicity of numeri- 
cal discretization and of handling boundary and interface conditions. A proper selection 
of the fundamental unknowns and a proper ordering of the governing equations can lead to 
further simplifications which produces a symmetric coefficient matrix in the governing 
equations. Of all the finite -difference schemes used with the first-order equation formu- 
lation, the scheme presented herein, which is based on the use of the two interlacing grids 
for the different fundamental unknowns, leads to the minimum bandwidth of the finite - 
difference field equations and the maximum accuracy of the solution and, therefore, is 
computationally most efficient. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., September 19, 1973. 
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APPENDIX A 


ELASTIC COEFFICIENTS OF LAMINATED SHELLS 
Elastic Stiffnesses of the Layers 

The stiffness coefficients of the kth orthotropic layer of the shell are given by: 


and 




E 


(k) 


;(k) 


Symmetric 


E (k) 
,(k) E L 

LT #) 

■p( k ) 
ill rp 


i(k) 

LT 


’i = 1,2,3\ 
j = l,2,3j 


(Al) 


[ c « 3 ■ 


j(k) 

TT 


G' 


(k) 

LT 


(!:«) 


(A2) 


where subscripts L and T denote the direction of fibers and the transverse direction, 
^ LT is the Poisson’s ratio measuring the strain in the T-direction due to a uniaxial 
normal stress in the L-direction, 


^LT E L " l/ TL E T 




(A3) 


x - 1 - ^ LT ^ TL 


J 


and superscript k refers to the kth layer. 


17 


APPENDIX A - Concluded 


Elastic Coefficients of the Shell 
The dimensionless elastic stiffnesses of the shell are given by: 



i] 




(A 4) 


and 


c ij 


C c k=l 



dx 3 



(A5) 


where 

NL = total number of layers of the shell 

h k and h k _i are the distances from the reference surface to the top and bottom 
surfaces of the kth layer, respectively 

C Q is a reference extensional rigidity of the shell 

The dimensionless elastic compliances of the shell Ay, By, Gy (i,j = 1,2,3) and 
Aij (i,j = 4,5) are obtained by inversion of the matrix of the elastic stiffnesses as follows: 


A ! B 

t 

i 


C ! F 
1 

i 

_B j G_ 


_F i D_ 


The elastic stiffnesses and compliances 
coordinate <f. 


M ■ ft ]' 1 

the shell are functions of the meridional 
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APPENDIX B 


FORMULAS FOR COEFFICIENTS IN GOVERNING 
DIFFERENTIAL EQUATIONS 

The independent nonzero terms of the submatrices [s a ] , [s b ] ? and [s] in equa- 
tions (17) are given by: 

b-, 


5 11 


P 


q a b 2 
b 12 ~ V 


5 13 


mb c b 3 


S 14 " ‘ ^ b c b 3 ' 


5 15 


mb c b 4 


3 22 


*23 


mb b e 
c 6 


3 24 


'e b c b 6 


3 25 


*33 


mb c b 7 

P 

2 2 
m b c G22 P^d ^2 


p 


□5 


S34 = -mwg | b c G 22 + i*. 



3 a 

5 35 


m2b c B 22 | P“e l 


v 55 
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APPENDIX B - Continued 


a 2 m k 2 

S 44 " b c^ G 22 p -- T 

^ A 55 

k 2 \ 


S7 C = m b B 


J 45 


c D 22 u e A 


55 


’55 


m b c A 22 P k 2 


55 


c b u 2 r 

S 11 - 'V 7 G 22 


S 12 “ 'V B 22 


c b _ m 
b i3 - 7 


o b 

S 14 ■ 7 


-m 


’15 


te ~ 7) 


S 22 _ " b c P7 A 22 


S b - -1 
b 24 " 1 


c b _ m 
b 25 - 7 


-b A 33 


33 


' \ {^6 ~^) A 33 + B 33 


35 
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APPENDIX B - Continued 



S H = V b 3 
S 12 = V b 4 

S 21 = V b 6 


5 31 = - m V G 22 

5 32 = m V B 22 
s 33 = y 

S35 = -y[^Q - <»g) 

5 41 = ' b c pyG 22 “d 

5 42 = b c pyB 22^ 

5 51 = m V B 22 

5 52 = " m V A 22 
S 55 = y 
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APPENDIX B — Continued 

where to are dimensionless coefficients given by: 


b l = A ll + b c 

{" A 12 ( A 12 G 22 “ B 12 B 22 

) + B 12 ( A 12 B 22 " A 22 B 12) 

b 2 = B 11 " b c 

{ A 12 ( B 12 G 22 " B 22 G 12) 

+ B 12 ( A 22 G 12 " B 12 B 22)} 

b 3 = A 12 G 22 " 

B 12 B 22 



b 4 = " A 12 B 22 + A 22 B 12 

b 5 = G 11 ' b c { B 12 ( B 12 G 22 ' B 22 G 12 ) + °12 ( A 22 G 12 " B 12 B 22)} 
b 6 = B 12 G 22 " B 22 G 12 
b 7 = “ B 12 B 22 + G 12 A 22 


with 


b c - ( A 22 G 22 “ B 22) 

and and ^ are tracing constants (shear coefficients). 
The matrix |l J is a 5 x 5 diagonal submatrix given by: 


W- 


-1 


The nonzero terms in the symmetric mass matrices [m a ] and [m b ] are 
given by : 

m 33 = k m 0 


nigg — k ^ k 
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APPENDIX B - Concluded 


and 


m 44 = k m 0 


m 55 = k r k m 2 


m 


n = k m o 


where 


m 12 = k r k m l 


m 22 = k r k m 2 


k=p^ 


k r is a tracing constant and 


NL h, 

L m ° m i m i\ L 

k=l k_1 

with p^ being the mass density of the kth layer of the shell 
s 

The components of the load vectors are given by: 
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TABLE 1 .- CHARACTERISTICS OF FINITE -DIFFERENCE SCHEMES 
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equations reduce to eight first-order equations. 



TABLE 2.- CONVERGENCE OF MINIMUM NONDIMENSIONAL FREQUENCY X 
OBTAINED BY DIFFERENT FINITE-DIFFERENCE SCHEMES 

jciamped isotropic cylindrical shell with h/R = 0.002, 

L/R = 10, v = 0.3, and m = 4; X - 100wR \J p s (\ - r 2 )/E 


n 

Values of x for - 

SW and DW 

MW 

MHS 1 and MHS 2 

MHI 

10 

2.101 

1.398 

1.549 

1.485 


a (1.393) 

(0.665) 

(1.027) 

(0.985) 

20 

1.689 

1.485 

1.517 

1.506 


(1.120) 

(0.985) 

(1.006) 

(0.999) 

30 

1.593 

1.501 

1.512 

1.508 


(1.056) 

(0.995) 

(1.003) 

(1.000) 


a Numbers in parentheses refer to 


converged- 
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TABLE 3.- ERROR INDEX FOR STRESS RESULTANTS AND 
GENERALIZED DISPLACEMENTS OBTAINED BY 
DIFFERENT FINITE-DIFFERENCE SCHEMES 









TABLE 4.- CONVERGENCE OF MINIMUM NONDIMENSIONAL FREQUENCY 
OBTAINED BY DIFFERENT FINITE-DIFFERENCE SCHEMES 


Clamped graphite-epoxy spherical shell with h/R = 0.002; m = 2; = cu ^h^/Er 




Values of 

x 10^ for - 


Converged 

n 

MW 

MHS 1 

MHS 2 

MHI 

solution 

6 = — (hemispherical shell) 
2 

10 

2.7962 

2.7830 

2.8614 

2.8403 

2.8217 

20 

2.8099 

2.7909 

2.8302 

2.8257 


30 

2.8146 

2.8095 

2.8266 

2.8239 


40 

2.8166 

2.8136 

2.8246 

2.8231 


| CO 
II 

Qi 

10 

3.4550 

— 

3.6456 

3.7724 

3.6489 

3.6644 

20 

3.6140 

3.6789 

3.6882 

3.6614 


30 

3.6411 

3.6704 

3.6760 

3.6634 


40 

3.6508 

3.6666 

3.6713 

3.6640 
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Figure 1.- Shell geometry. 




Figure 2.- Sign convention for stress resultants and displacements. 
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Figure 3.- Control points for the different fundamental unknowns in the MHI scheme. 




(a) MW scheme. 

Figure 4.- Schematic representation of the finite -difference equations 
obtained by the first-order equation schemes. 











(b) MHI scheme. 
Figure 4.- Continued, 
































h 



h/R * 0.002 
L/R - 10 
v * 0.3 

Figure 5.- Characteristics of isotropic cylindrical 
shell used in the present study. 
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Meridional distance, £ 


(a) Normal displacement, W. 

Figure 6.- Mode shapes and modal stress resultants obtained by MHS 1, MHS 2 
and MHI schemes for isotropic circular cylindrical shell. 







ratio 








P r * p 0 cos m 6 


Figure 7.- Characteristics of toroidal shell and loading. 



Meridional stress resultant, N^/(p o R) Normal displacement, WE/(p h) 










Number of finite-difference intervals, n 


Figure 9.- Convergence of minimum frequency obtained by different finite- 
difference schemes based on the first-order equation formulation 
(graphite-epoxy toroidal shell with h/R = 0.02; m = 2). 
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Fibers of top layer in the meridional direction 




Subscripts L and T refer to direction of 
fibers and transverse direction, 
respectively. 

Figure 10.- Eight-layered graphite -epoxy spherical shell with clamped edges. 
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